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We use a tight-binding Bogoliubov-de Gennes (BdG) formalism to self-consistently calculate the 
proximity effect, Josephson current, and local density of states in ballistic graphene SNS Josephson 
junctions. Both short and long junctions, with respect to the superconducting coherence length, 
are considered, as well as different doping levels of the graphene. We show that self-consistency 
does not notably change the current-phase relationship derived earlier for short junctions using the 
non-selfconsistent Dirac-BdG formalism^ but predict a significantly increased critical current with 
a stronger junction length dependence. In addition, we show that in junctions with no Fermi level 
mismatch between the N and S regions superconductivity persists even in the longest junctions we 
can investigate, indicating a diverging Ginzburg-Landau superconducting coherence length in the 
normal region. 

PACS numbers: 74.45.+C, 74.50.+r, 74.78.Na, 73.20.At 



I. INTRODUCTION 

Recently SNS graphene Josephson junctions have been 
experimentally realized by depositing s-wave supercon- 
ducting contacts on top of a graphene sheet^^ These 
junctions have been shown to carry a Josephson current 
that, depending on the position of the Fermi level relative 
to the Dirac point, consists of either electrons or holes, 
and, even at the Dirac point, where the density of states 
is zero, a finite supercurrent has been measured. 

Theoretically, a finite Josephson current was predicted 
by Titov et. a& shortly before the first experimental real- 
ization. This result was based on the Dirac-Bogoliubov- 
de Gennes (DBdG) formalism developed by Beenakker— 
where the band structure is approximated by the par- 
ticular double Dirac cone spectra found in graphene and 
the traditional BdG formalism is then applied to solve for 
the wave functions throughout the junction. The Joseph- 
son current was approximated as the current carried by 
the subgap bound Andreev states in the junction, an ap- 
proximation valid in the short junction regime^ where 
the length L of the junction is smaller than the super- 
conducting coherence length £. Graphene SBS junctions, 
where B is a barrier created by a heavily doped graphene 
slip, were later also investigated using the same general 
formalism.— 

The DBdG approach in Ref. Q neglects the spacial de- 
pendence of the superconducting order parameter as it 
assumes a constant, non-zero, order parameter in the su- 
perconductor and an abrupt change to zero at the SN 
interface. However, close to a SN interface the supercon- 
ducting order parameter is expected to vary strongly as 
a function of the distance to the interface. In this paper 
we address this problem by self-consistently calculating 
the superconducting order parameter in graphene SNS 
Josephson junctions using the tight-binding Bogoliubov- 



de Gennes (TB BdG) formalism. Specifically, this allows 
us to explicitly calculate the proximity effect depletion 
of the order parameter close to the interface and it also 
results in a Josephson current properly calculated from 
the superconducting proximity effect. In addition, the 
results will no longer be limited to the short junction 
regime. 

As in the DBdG approach, we will study an impurity 
free, ballistic, graphene sheet by itself, and therefore need 
to model the influence of the superconducting contacts on 
the graphene. We do this by assuming the following two 
effects on the graphene directly under the superconduct- 
ing contacts: The superconducting contacts induce 1) 
an on-site attractive Hubbard pairing potential U, which 
will lead to a s-wave superconducting order parameter 
Ay, and 2) a heavy doping in the graphene. The first as- 
sumption implies that the superconducting contacts in- 
duce a pairing potential in the graphene but that the 
superconducting order parameter itself in the graphene 
junction will be subjected to the proximity effect when 
solved for self-consistently. Since solving self-consistently 
not only for the order parameter but also for the chem- 
ical potential in the system adds a significant amount 
of complexity to the problem, we will assume that the 
doping profile throughout the junction is set by differ- 
ent, but constant, effective chemical potentials fx in the 
S and N regions. The same assumption was used in the 
non-self-consistent DBdG approach. Figure [1] shows a 
schematic of the experimental (a) and model setup (b) . 
The interface between S and N is assumed to be clean 
and smooth, a realistic assumption as experiments has 
indicated a high transparency of the SN interfaces^. 

Using our self-consistent TB BdG formalism we find 
that the functional dependence of the supercurrent on 
the phase drop across the junction is not significantly al- 
tered from the non-self-consistent DBdG results by Titov 
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FIG. 1: Schematic of (a) an experimental SNS graphene 
Josephson junction and (b) the model setup with input pa- 
rameters: pairing potential U, effective chemical potential ft, 
length of normal region L, and region where the phase of the 
order parameter will be kept fixed X. 



et. aJi and in fact extends to junctions with L > £. How- 
ever, we find that the absolute value of the critical cur- 
rent is in general significantly higher and has a stronger 
L-dependence, especially for junctions where the N re- 
gion is moderately doped. In addition, we show that 
for junctions with no Fermi level mismatch (FLM) be- 
tween the N and S regions superconductivity persists, 
not just proximity- wise, even in the longest junctions we 
can study. Finally, we also report detailed local den- 
sity of states (LDOS) plots for different SN junctions, 
showing the evolution of the superconducting energy gap 
throughout the junction. These results should directly 
be comparable with experimental data obtained using a 
point contact scanning tunneling probe. 

Before proceeding it is worth noting that the TB BdG 
formalism not only effectively captures the proximity ef- 
fect in a s-wave graphene Josephson junction but is also 
easily extendable to include other short-range electronic 
coupling terms in the graphene. Going back originally 
to Linus Pauling^ jwr-bonded planar organic molecules, 
such as graphene, have been recognized to have enhanced 
nearest-neighbor spin-singlet bonds compared to polar 
configurations. In a Hamiltonian formulation this spin- 
singlet enhancement takes the form of an intrinsic JSi ■ Sj 
term between nearest neighbors which for strong enough 
coupling will give rise to mean-field superconductivity. In 
earlier work^ we have shown that in the bulk this gives 
rise to a time- reversal symmetry breaking d x 2_ y 2 + id xy 
superconducting order parameter. In addition, a JSi ■ Sj 
term can also be used to model (i-wave superconducting 
contacts in graphene Josephson junctions. Self-consistent 



studies including these intrinsic spin-singlet electronic 
correlations in the graphene inside a s-wave SNS Joseph- 
son junction as well as the effect of d-wave contacts are 
the subject of future publications. 



II. METHOD 

Based on the motivation above we model a graphene 
SNS Josephson junction using the following effective 
tight-binding, attractive Hubbard Hamiltonian 

<i,j>,er i.er 

- E CWj/iT/i/u + <h<>v<h.<n/- (1) 



Here fy is the creation operator on the A-site in cell 
i = (n, to) of the honeycomb lattice, and g\ on the B-site 
of the same unit cell i, see Fig. [2] <i, j > indicates a sum 
over nearest neighbors. The energy parameters are the 
hopping energy t w 2.5 eV which we will assume constant 
throughout the junction, the on-site attractive Hubbard 
term U (i) which is only non-zero in the S regions, and the 
effective chemical potential which includes the com- 
bined effect of the chemical potential /i, which is constant 
in the whole system, and the local doping. To proceed, 




(n+l,m-3) 



FIG. 2: (Color online) The graphene honeycomb lattice with 
the two different atomic sites / and g, the unit cell vectors 
{ci,C2}, the three nearest neighbor directions {ai,a2,a3}, 
the zigzag interface with its T = a = 2.46 A long unit cell 
in the y-direction, and the (n, m) notation for labeling each 
unit cell. 

we use the Hartree-Fock-Bogoliubov mean-field approx- 
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imation to arrive at the following one-particle Hamilto- 
nian 

H MF = - *E Ulsi,* + slfi,*) +E A(i)(/iU* + glgia) 

<i,j>,cr i.a 

+ EA [/ (i)(// T / i t i+ .gt T5 t i)+H .c., (2) 



where we have defined the spatially dependent mean-field 
superconducting order parameter 



Au(i) = -U(i) {fM+ J 9im) . 



(3) 



It is straightforward to show that the standard tight- 
binding Bogoliubov-de Gennes (BdG) formalism (see e.g. 
Refs. llOUllUl2l for recent applications) applies even in the 
case of the 2-atom unit cell found in graphene. Writ- 
ten out explicitly, the mean-field Hamiltonian |(5J) can 
be diagonalized with the following site-dependent two- 
dimensional Bogoliubov-Valatin transformation 



2iV 




(4) 



(5) 



where TV is the number of unit cells in the whole junction, 
and the resulting diagonal Hamiltonian reads 



2N 



Hmf 



(6) 



The energy eigenvalues E v and the eigenfunctions 
(u, y, v, z)\ are determined by solving the following 47V x 
47V eigenvalue problem for its 27V positive eigenvalues 



E 



Ho(i,j) A(i,j) 
At(iJ) -ff (i,j) 



/ u l \ 
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(7) 



where Hq and A are 2x2 matrices which for (i,j) 
((n, m), (p, r)) can be written as 



#o(i,j) = 



A(i)<5ij 



-iS, 



r,m=pl / 



/i(i)$y 



. ( A^(i)J u 



A(U) = 











(8) 
(9) 



As in the standard 1-atom per unit cell BdG formalism, 
half of the eigenvalues are guaranteed to be positive be- 
cause Hq = Ho and A = A T . The self-consistency con- 

>e writ 

f5E v 



dition for the order parameter can finally be written as 
Mi) = +^Dtanh^. ( ID) 



To significantly reduce the size of eigenvalue problem, 
we will assume a smooth interface between the S and N 
regions and use the translational symmetry that exists 
perpendicular to the junction to apply Bloch's theorem. 
The formalism and subsequent numerical results are for 
the zigzag interface (see Fig. (2) but we have also studied 
the armchair interface and found no differences. With 
only on-site superconducting pairing the direction of the 
interface should indeed not matter. More specifically, for 
the zigzag interface we can now write the eigenfunctions 
as 



f w i— (n.m) 
V\— (n,m) 
^i— (n,m) 

V Z'l— (n,m) 



E 



Un{ky) 
V z n(ky) 



(11) 



where the wave vector k v 



~ TTa' w ith I being an integer 
such that k y s] — , — ], is a good quantum number for the 
system and 7V y denotes the number of unit cells, of width 
T = a = 2.46 A, perpendicular to the interface. The BdG 
eigenvalue problem Eq. ([7|) now reduces to only depend 
on the x coordinate indices (n,p) but has to be solved 
for all k y . 

The BdG eigenvalue equations above, Eq. ((7l l9ll0p are 
solved by starting with an initial guess for the order pa- 
rameter profile Ajj throughout the junction. After find- 
ing the 27Va; eigenstates of Eq. ([7]) for each k y we can 
compute a new order parameter profile from Eq. (1 10[) . 
The process is repeated until the difference in order pa- 
rameters between two subsequent steps are less than a 
desired accuracy. The final converged result will allow 
us to study the proximity effect in graphene. While Ajj 
will always be zero outside S, since per definition it is 
proportional to U, the pairing amplitude 



FuO) = 



(/u/h) + (sum) 



Mi) 

U(i) 



(12) 



will display the leakage of Cooper pairs into the normal 
region. 

In the above formalism it is the leakage pairing ampli- 
tude that is responsible for the Josephson current when 
the order parameter has a finite phase gradient over the 
junction. Numerically, we impose a phase gradient by 
fixing the phase of Ajj in the outermost parts of the S re- 
gions, labeled X in Fig.[TJ and then solve self-consistently 
for both phase and amplitude in the remaining SNS junc- 
tion but only for the amplitude in the X regions. For 
small currents one should be able to fix the phase in the 
whole S regions, but physically, as soon as the current 
is non-zero, there will necessarily also be a phase drop 
even over the superconducting contacts and not just in 
the normal region. By extensive testing we have found 
this to be a non-insignificant effect in many junctions 
and we have therefore ensured that the region of self- 
consistency for the phase in the contacts is large enough 
to ensure bulk-like conditions in S. In terms of calculat- 
ing the Josephson current vs. phase drop, 1 ((/>), we still 
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define the phase variable <fi as the phase drop over the 
normal region N as is usually done. 

Once a self-consistent solution with a phase drop <p 
over the junction is obtained, we calculate the Joseph- 
son current using the continuity equation for the particle 
current 

V.J + |=0 (13) 

together with the Heisenberg equation 

dm i , , 
-^ = - h \HM, d4) 

where m is the particle density per unit cell with the av- 
erage p = (n). This approach was used in e.g. Ref. [H for 
a square lattice. The quantum average of the commuta- 
tor in Eq. (fl4|) can easily be shown to only contain the 
kinetic hopping terms when the self-consistent order pa- 
rameter is used in the mean- field Hamiltonian @ , which 
is true in our model except in the X regions. The X 
regions will therefore act as sinks and/or sources for the 
current. More specifically, we get the Josephson current 
per cross-sectional distance as 

I(n) = ea cc i\J (n) / (Ha KS ) (15) 

where 

Of 

J » = " W E M< t (*y)l/n-i(A»))/(^") 

+ Im«(fe y )<L 1 (fc 2/ ))/(-^)) cos(fcW2) 

(16) 

and, for the zigzag interface, a ce u = \/3a/2 is the length 
of the unit cell in the direction of the current and a xs = 
a/2 is the perpendicular, cross-sectional, distance. 

Finally, we will also be interested in the local density of 
states (LDOS) which easily can be calculated at T = K 
from 

D i (E) = 2j2(\Ui\ 2 + \yn 2 )S(E"-E) 

+ (M\ 2 + W\ 2 )S(E + En. (17) 

III. NUMERICAL RESULTS 

A. Simulation details 

To investigate the proximity effect and Josephson cur- 
rent in graphcnc SNS junctions we have solved the tight- 
binding BdG formalism described above self-consistcntly. 
The physical input parameters arc the on-site pairing po- 
tential U in S, the effective potential p, in S and N, tem- 
perature, and the length L of N. More specifically we 
consider the following set up. For the superconducting 
contacts: U(S) = 3.4 eV = 1.36t, fi(S) = 1.5 eV = 0.6t. 



This leads to Ay = Ao = 0.1 eV and a superconducting 
coherence length £ « 50 A ss 25 unit cells in the zigzag 
direction for the bulk. These values satisfy Xf(S) <C £ 
and allow us to numerically investigate both the L < £ 
and L > £ cases. They are however are quite large val- 
ues for a realistic situation but smaller superconducting 
gaps leads to slower convergence rates and also the need 
for larger systems making calculations less feasible. We 
have checked our key results for smaller U and found 
no significant difference. For the normal region we have 
mainly studied short junctions with L = 10 unit cells and 
long junctions with L = 50 unit cells with various doping 
levels. The doping levels have been implemented by set- 
ting n(N) to values ranging from eV, modeling clean, 
undoped, graphene where the Fermi level is located at 
the Dirac point to heavily doped junctions with 1.5 eV 
which gives no FLM at the interfaces. The temperature 
was chosen to be T — 10 K throughout the work, which 
in comparison with T c is effectively K. 

The accuracy of the solution is determined by the 
choice of termination criterion for the self-consistency 
step, the number of k y points, and the size of S and X to 
ensure bulk- like superconducting conditions, all of which 
has been tested thoroughly. 



B. Proximity effect 

Fig. [3] shows typical, normalized, pair amplitude Fjj 
profiles through a SNS junction with different doping 
levels in N. The depletion of superconducting pairs in 
S close to the junction as well as the leakage of pairs into 
the junction is clearly visible. The oscillations at each 
end of S are due to the end surfaces of the contacts, and 
are not of primarily importance here. For small doping 
levels, there are also pronounced oscillations in Fjj at 
the interface on the S side. These oscillations are cor- 
related with oscillations in the charge density which are 
present even for the case U = and therefore attributed 
to Friedel-like charge oscillations due to the FLM at the 
interface. 



C. I vs. curves 

As Fig. [3] clearly shows there exists a pronounced prox- 
imity effect in graphene SNS junctions. This effect is 
obviously not taken into account in a calculation where 
the order parameter is assumed to change abruptly at 
the interface, such as in Ref. [l|, and the natural ques- 
tion arises if this effect will significantly change any pre- 
dictions made by such a non-self-consistent calculation. 
Among the most significant quantities of a SNS junc- 
tion is the Josephson supercurrent it can sustain if the 
junction is short enough to allow coherent transport of 
superconducting pairs. We investigate in the following 
two subsections two properties of the Josephson current. 
In this subsection we extract the current vs. phase re- 




L= 10 



L = 50 



Unit cell 

FIG. 3: (Color online) Normalized pair amplitude Fu profiles 
for a S = 50, N = 50 unit cells junction with three different 
doping levels in N: at Dirac point eV (black), moderate 
doping 0.7 eV (red), no FLM 1.5 eV (green). 



lationship I{4>) for junctions with different lengths and 
doping levels and in the next subsection we will look at 
the length dependence of the critical current for short 
junctions. Both of these results were extensively worked 
out in the DBdG non-self-consistent formalism in Ref. 1. 

Before proceeding it should be noted that the largest 
phase drop we can put over the whole structure is tt. 
However, with a non-zero phase drop over the S regions 
the maximum phase drop <j) over the junction will be 
smaller than tt for junctions with a significant current and 
we will in these cases not be able to numerically trace out 
the I(4>) relationship over the full [0, tt] interval. While 
this appears as a numerical artifact in this context it is in 
fact closely related to the physical 2tt phase-slip process 
in Josephson junctions (see e.g. Ref. fl3l). 

Fig. |4] shows the phase dependence of the Josephson 
current for both short junctions (L < £) and long junc- 
tions [L > £) for three different doping levels; undoped 
p.{N) — eV, moderately doped such that p.(N) = 
0.7 eV ^> Hvf/L, and with no FLM at the interface. The 
currents have been normalized to the maximum value of 
the undoped, short junction. The self-consistent data 
are in black whereas the red and green curves are least 
squares fits to results obtained in the DBdG formalism. 
More specifically, for the undoped junction the red curve 
is a fit, with C as the fitting parameter, to Eq. (20) in 
Ref. [S 



I(4>) = Ccos(<5y2)arctanh[sin(</>/2)]. 



(18) 



This phase relationship is notably different from the tra- 
ditional Josephson formM, I — I c sin(0) , but is in fact 
identical to that of a disordered metal upon substitution 
kpl — > 1, where / is the mean free path. This is despite 
the fact that the graphene SNS junction is treated as a 
ballistic junction and is here instead a consequence of 




FIG. 4: (Color online) I(<f>) relationsship, normalized with re- 
spect to the critical current in (e) for short L < £ (a,c,e) and 
long L > £ (b,d,f) junctions when the N region is doped at 
the Dirac point jl — eV (e,f), moderately doped jl = 0.7 eV 
(c,d), and with no FLM jl = 1.5 eV (a,b). Self-consistent 
numerical data are in black with the lines being guides to 
the eye. Red curves are least squares fits (of the overall con- 
stant) to the DBdG results whereas green curves are fits to 
the functional form sgn(cos(</>/2)) sin(^/2). 



the Dirac band spectrum of graphene. As can be seen 
in Fig. [4] there is good agreement to this functional de- 
pendence for the short junction at the Dirac point (e), 
but also for the long junction (f) where the approxima- 
tion of calculating the Josephson current from the sub- 
gap Andreev bound states in the DBdG formalism is not 
formally motivated. For non-zero, but moderate, dop- 
ing of the junction a closed form analytical result from 
the DBdG formalism is not available, though the func- 
tional dependence on cf> is still close to Eq. [T8l and there 
is still good agreement with our self-consistent results 
(Fig. Q]c,d). Finally, for the case of no FLM between 
the S and N regions (Fig. [5] a, b), the green curves are fit 
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to the functional form sgn(cos(</>/2)) sin(0/2) which was 
derived recently by Linder et. afe Due to the high cur- 
rents in the no FLM cases, it is not possible to obtain cur- 
rent values for phase drops larger than those reported in 
Fig. [4l and it is thus impossible to distinguish which fits 
are best. In fact, for both of these two junctions as well 
as for the longest junctions we could conveniently model 
(L 210 A) with no FLM, the superconductivity energy 
gap is never depleted inside the junction, and, subse- 
quently, the phase drop profile across the whole structure 
is linear. A linear phase drop across the whole structure 
will give a linear relationship between / and </> which is 
clearly seen in Fig. |4j This persistent superconductivity 
gap in the junction has the intriguing consequence that 
for these junctions the normal region has, at least in our 
structure sizes, a diverging Ginzburg-Landau coherence 
length, despite the fact that the pairing potential is zero 
inside the junction. 



D. L dependence of the critical current 

In the previous subsection the I vs. <j> dependence was 
investigated and fitted to analytical forms derived from 
the non-self-consistent DBdG formalism. No attention 
was however paid to the dependence of the prefactor, 
i.e. the fitting parameter, or equivalently the critical cur- 
rent, on various physical quantities such as doping level 
in N, junction length, gap size etc. In this subsection we 
will focus on the behavior of the critical current I c on 
the length of the junction L in the short junction regime. 
More specifically, we will investigate the influence of the 
proximity effect on the L dependence of the following two 
results derived using the DBdG formalism 1 



eA 1 

I c = 1.33— when jl = 0, 

n ttL 



I c = 1.22 



eA fi 
H ttHvf 



(19) 

when jl > hvp/L. (20) 



Fig. E] shows the the I c vs. L dependence for the self- 
consistent results (crosses) and Eqs. (|19H20[) (black), to- 
gether with least squares fits to the functional forms CL b 
(red) and Ce~ L ^ N ( green). The current is normalized 
such that I C (L = 10) = 1 for Eq. (f2"0"|) . For a junction 
at the Dirac point, Fig. [5]Ja), we see that the proximity 
effect causes an increase in the critical current as well 
as change the L dependence from b = — 1 in Eq. (fT9j) to 
b = —1.3. Functionally this form is also close to the tradi- 
tional Ginzburg-Landau Ce~ L ^ N functional dependence 
with £jv = 16 A. At moderate doping, here represented 
by jl = 0.7 eV 3> Hvf/L, Fig.[5jb) shows that the effect 
of a self-consistent calculation is even more pronounced. 
Eq. (f2"0"|) has no L dependence whereas we see a clear 
increase in the critical current when the junction size de- 
creases. The fits are not as good as in the undoped case 
but b « —0.4 or £jy ~ 50 A. Also note the significant in- 
crease, between 2 and 4 times, in the current for all junc- 
tion lengths investigated. Interestingly, the difference in 




5 10 15 
L (unit cells) 



10 15 
L (unit cells) 



FIG. 5: (Color online) Normalized critical Josephson current 
I c as function of length L of the junction at (a) the Dirac point 
jl = eV and (b) moderately heavy doping ji — 0.7 eV. Self- 
consistent results (crosses), Eqs. (|19I2U|I (black), and least- 
square fits to the functional forms CL b (red) and Ce~ L ^ N 
(green). The current is normalized such that the critical cur- 
rent in Eqs. (|20p is equal to 1. 



6-exponents between the two cases are ~ 1, just as the 
DBdG prediction, although both values are significantly 
modified in our self-consistent calculations. 



E. LDOS for SN junction 

Finally, we report in Fig. [5] on how the local density of 
states (LDOS) evolve through a graphene SN interface 
with an abrupt effective chemical potential jl change at 
the interface. Since for a non-uniform system, such as 
a SN junction, the superconducting energy gap will in 
general not be equivalent to the order parameter Ay in 
the full self-consistent solution we are especially inter- 
ested in the evolution of the superconducting energy gap 
and the adjoining coherence peaks. The relevant super- 
conducting length scale is the coherence length which is 
£ « 50 A w 25 unit cells on the S side. The two up- 
per plots (a,b) show the LDOS when the normal region 
is undoped, with the right hand plot being a magnifica- 
tion around zero energy, defined to be at the Dirac point. 
As can be seen, the superconducting coherence peaks on 
each side of the energy gap are pronounced, and in fact 
enhanced close to the interface, all the way up to the in- 
terface but die out very quickly, over only a few unit cells, 
on the N side. Since the DOS is zero at the Dirac point it 
is hard with the current resolution to determine the ex- 
act evolution and completeness of the superconducting 
energy gap. For the two lower plots (c,d) the N region 
is instead biased into a moderately heavy doping regime. 
Here the superconducting gap is not complete even at 
distances ~ 40 A from the interface on the S side. Close 
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E/A () E/A 

FIG. 6: (Color online) Typical LDOS for SN junction as func- 
tion of the normalized energy denned with reference to the 
Dirac point. Upper plots (a,b): N region is undoped, at the 
Dirac point. Lower plots (c,d): N is moderately doped. Right 
hand plots are magnifications of the gap features of the left 
hand plots. Different colors represent different positions in- 
side the junction: middle of S (black), middle of N (red), 
interface (green), interface ±2 unit cells (magenta, cyan), in- 
terface —20 unit cells (blue). 

to the interface the coherence peaks are shifted towards 
lower energies. On the N side of the interface a normal 
state LDOS is again achieved within only a few unite 
cells. 

The above results thus give a very short decay length 
for the superconducting state in the N region but show 
on large and spatially extended effects of the supercon- 
ducting gap on the superconducting S side, especially as 
we move away from the Dirac point. This should be con- 
trasted with the results for junctions with no FLM where 



a superconducting state, with an complete energy gap, 
persist in N even in the longest junctions we can model. 
The significant different response between undoped or 
moderately doped junctions on one hand and junctions 
with no FLM on the other hand should be experimen- 
tally accessible using a point contact scanning probe to 
investigate the LDOS. 

IV. CONCLUSIONS 

Using the TB BdG formalism we have been able to cal- 
culate the proximity effect and Josephson current in bal- 
listic SNS graphene Josephson junctions. We have shown 
that the functional form of /(</>) derived using the non- 
self-consistent DBdG formalism for short junctions^ is 
qualitatively valid at both in undoped and in moderately 
doped junctions. In addition, we have demonstrated that 
these results are extendable to the long junction regime 
where the junction length L is longer than the super- 
conducting coherence length. However, the dependence 
on the junction length L for the critical current in short 
junctions is enhanced, from L^ 1 to L~ 13 for undoped 
and from L° to L~ ' 4 for moderately doped junctions, 
when proximity effect is taken into account. Also, the 
magnitude of the critical current is enhanced. For junc- 
tions with no FLM between S and N, i.e. the same doping 
level throughout the structure, superconductivity is not 
depleted inside N even for the longest junctions we were 
able to investigate. This means that in a ballistic, no 
FLM, graphene SNS junction the Ginzburg-Landau su- 
perconducting coherence length in the normal region is 
close to diverging. This intriguingly indicates that the 
normal region is close to a superconducting instability, 
despite the fact that no pairing potential is present in 
this region. This behavior is in sharp contrast with the 
LDOS data for SN junctions when N is undoped or only 
moderately doped where the superconducting gap disap- 
pears within only a few unit cells on the N side of the 
junction. 
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